// opengl visualization/gui for code found on Paul Bourkes page:
// http://local.wasp.uwa.edu.au/~pbourke/papers/triangulate/
// its the C++ version by Gilles Dumoulin
#pragma comment(lib,"opengl32.lib")
#pragma comment(lib,"glfwdll.lib")

#include "Delaunay.h"
#include <gl\glfw.h>
#include <vector>
#include <sstream>

#include "misc_utils.h"


using namespace std; 

////////////////////////////////////////////////////////////////////////
// CircumCircle() :
//   Return true if a point (xp,yp) is inside the circumcircle made up
//   of the points (x1,y1), (x2,y2), (x3,y3)
//   The circumcircle centre is returned in (xc,yc) and the radius r
//   Note : A point on the edge is inside the circumcircle
////////////////////////////////////////////////////////////////////////

int CircumCircle(double xp, double yp, double x1, double y1, double x2, 
double y2, double x3, double y3, double &xc, double &yc, double &r){
  double m1, m2, mx1, mx2, my1, my2;
  double dx, dy, rsqr, drsqr;

/* Check for coincident points */
if(abs(y1 - y2) < EPSILON && abs(y2 - y3) < EPSILON)
  return(false);
if(abs(y2-y1) < EPSILON){ 
  m2 = - (x3 - x2) / (y3 - y2);
  mx2 = (x2 + x3) / 2.0;
  my2 = (y2 + y3) / 2.0;
  xc = (x2 + x1) / 2.0;
  yc = m2 * (xc - mx2) + my2;
}else if(abs(y3 - y2) < EPSILON){ 
        m1 = - (x2 - x1) / (y2 - y1);
        mx1 = (x1 + x2) / 2.0;
        my1 = (y1 + y2) / 2.0;
        xc = (x3 + x2) / 2.0;
        yc = m1 * (xc - mx1) + my1;
      }else{
         m1 = - (x2 - x1) / (y2 - y1); 
         m2 = - (x3 - x2) / (y3 - y2); 
         mx1 = (x1 + x2) / 2.0; 
         mx2 = (x2 + x3) / 2.0;
         my1 = (y1 + y2) / 2.0;
         my2 = (y2 + y3) / 2.0;
         xc = (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2); 
         yc = m1 * (xc - mx1) + my1; 
       }
  dx = x2 - xc;
  dy = y2 - yc;
  rsqr = dx * dx + dy * dy;
  r = sqrt(rsqr); 
  dx = xp - xc;
  dy = yp - yc;
  drsqr = dx * dx + dy * dy;
  return((drsqr <= rsqr) ? true : false);
}
///////////////////////////////////////////////////////////////////////////////
// Triangulate() :
//   Triangulation subroutine
//   Takes as input NV vertices in array pxyz
//   Returned is a list of ntri triangular faces in the array v
//   These triangles are arranged in a consistent clockwise order.
//   The triangle array 'v' should be malloced to 3 * nv
//   The vertex array pxyz must be big enough to hold 3 more points
//   The vertex array must be sorted in increasing x values say
//
//   qsort(p,nv,sizeof(XYZ),XYZCompare);
///////////////////////////////////////////////////////////////////////////////

int Triangulate(int nv, XYZ pxyz[], ITRIANGLE v[], int &ntri){
  int *complete = NULL;
  IEDGE *edges = NULL; 
  //std::vector<*IEDGE> edges;
  IEDGE *p_EdgeTemp;
  int nedge = 0;
  int trimax, emax = 200;
  int status = 0;
  int inside;
  int i, j, k;
  double xp, yp, x1, y1, x2, y2, x3, y3, xc, yc, r;
  double xmin, xmax, ymin, ymax, xmid, ymid;
  double dx, dy, dmax; 

/* Allocate memory for the completeness list, flag for each triangle */
  trimax = 4 * nv;
  complete = new int[trimax];
/* Allocate memory for the edge list */
  edges = new IEDGE[emax];
/*
      Find the maximum and minimum vertex bounds.
      This is to allow calculation of the bounding triangle
*/
  xmin = pxyz[0].x;
  ymin = pxyz[0].y;
  xmax = xmin;
  ymax = ymin;
  for(i = 1; i < nv; i++){
    if (pxyz[i].x < xmin) xmin = pxyz[i].x;
    if (pxyz[i].x > xmax) xmax = pxyz[i].x;
    if (pxyz[i].y < ymin) ymin = pxyz[i].y;
    if (pxyz[i].y > ymax) ymax = pxyz[i].y;
  }
  dx = xmax - xmin;
  dy = ymax - ymin;
  dmax = (dx > dy) ? dx : dy;
  xmid = (xmax + xmin) / 2.0;
  ymid = (ymax + ymin) / 2.0;
/*
   Set up the supertriangle
   his is a triangle which encompasses all the sample points.
   The supertriangle coordinates are added to the end of the
   vertex list. The supertriangle is the first triangle in
   the triangle list.
*/
  pxyz[nv+0].x = xmid - 20 * dmax;
  pxyz[nv+0].y = ymid - dmax;
  pxyz[nv+1].x = xmid;
  pxyz[nv+1].y = ymid + 20 * dmax;
  pxyz[nv+2].x = xmid + 20 * dmax;
  pxyz[nv+2].y = ymid - dmax;
  v[0].p1 = nv;
  v[0].p2 = nv+1;
  v[0].p3 = nv+2;
  complete[0] = false;
  ntri = 1;
/*
   Include each point one at a time into the existing mesh
*/
  for(i = 0; i < nv; i++){
    xp = pxyz[i].x;
    yp = pxyz[i].y;
    nedge = 0;
/*
     Set up the edge buffer.
     If the point (xp,yp) lies inside the circumcircle then the
     three edges of that triangle are added to the edge buffer
     and that triangle is removed.
*/
  for(j = 0; j < ntri; j++){
    if(complete[j])
      continue;
    x1 = pxyz[v[j].p1].x;
    y1 = pxyz[v[j].p1].y;
    x2 = pxyz[v[j].p2].x;
    y2 = pxyz[v[j].p2].y;
    x3 = pxyz[v[j].p3].x;
    y3 = pxyz[v[j].p3].y;
    inside = CircumCircle(xp, yp, x1, y1, x2, y2, x3, y3, xc, yc, r);
    if (xc + r < xp)
// Suggested
// if (xc + r + EPSILON < xp)
      complete[j] = true;
    if(inside){
/* Check that we haven't exceeded the edge list size */
      if(nedge + 3 >= emax){
        emax += 100;
        p_EdgeTemp = new IEDGE[emax];
        for (int i = 0; i < nedge; i++) { // Fix by John Bowman
          p_EdgeTemp[i] = edges[i];   
        }
        delete []edges;
        edges = p_EdgeTemp;
      }
      edges[nedge+0].p1 = v[j].p1;
      edges[nedge+0].p2 = v[j].p2;
      edges[nedge+1].p1 = v[j].p2;
      edges[nedge+1].p2 = v[j].p3;
      edges[nedge+2].p1 = v[j].p3;
      edges[nedge+2].p2 = v[j].p1;
      nedge += 3;
      v[j] = v[ntri-1];
      complete[j] = complete[ntri-1];
      ntri--;
      j--;
    }
  }
/*
  Tag multiple edges
  Note: if all triangles are specified anticlockwise then all
  interior edges are opposite pointing in direction.
*/
  for(j = 0; j < nedge - 1; j++){
    for(k = j + 1; k < nedge; k++){
      if((edges[j].p1 == edges[k].p2) && (edges[j].p2 == edges[k].p1)){
        edges[j].p1 = -1;
        edges[j].p2 = -1;
        edges[k].p1 = -1;
        edges[k].p2 = -1;
      }
       /* Shouldn't need the following, see note above */
      if((edges[j].p1 == edges[k].p1) && (edges[j].p2 == edges[k].p2)){
        edges[j].p1 = -1;
        edges[j].p2 = -1;
        edges[k].p1 = -1;
        edges[k].p2 = -1;
      }
    }
  }
/*
     Form new triangles for the current point
     Skipping over any tagged edges.
     All edges are arranged in clockwise order.
*/
  for(j = 0; j < nedge; j++) {
    if(edges[j].p1 < 0 || edges[j].p2 < 0)
      continue;
    v[ntri].p1 = edges[j].p1;
    v[ntri].p2 = edges[j].p2;
    v[ntri].p3 = i;
    complete[ntri] = false;
    ntri++;
  }
}
/*
      Remove triangles with supertriangle vertices
      These are triangles which have a vertex number greater than nv
*/
  for(i = 0; i < ntri; i++) {
    if(v[i].p1 >= nv || v[i].p2 >= nv || v[i].p3 >= nv) {
      v[i] = v[ntri-1];
      ntri--;
      i--;
    }
  }
  delete[] edges;
  delete[] complete;
  return 0;
} 

void randomize(){
  srand( (unsigned int) (time_t) time(NULL));  
}


int XYZCompare(const void *v1, const void *v2){
  XYZ *p1, *p2;
    
  p1 = (XYZ*)v1;
  p2 = (XYZ*)v2;
  if(p1->x < p2->x)
    return(-1);
  else if(p1->x > p2->x)
         return(1);
       else
         return(0);
}

void outputtriangle(size_t &nv, XYZ p[], ITRIANGLE v[], int &ntri){
	glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
	glBegin(GL_TRIANGLES);
	for(int i = 0; i < ntri; i++)
	{
		glVertex2d( p[v[i].p1].x, p[v[i].p1].y );
		glVertex2d( p[v[i].p2].x, p[v[i].p2].y );
		glVertex2d( p[v[i].p3].x, p[v[i].p3].y );
	}
	glEnd();
}

int main(){
	
	if( glfwInit() != GL_TRUE)
	{
		printf("glfw failed\n");
		return 1;
	}

	int xres=500;
	int yres=500;
	glfwOpenWindow(xres,yres,0,0,0,0,0,0,GLFW_WINDOW);
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	glOrtho(0,xres,yres,0,-1,1);
	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();

  ITRIANGLE *v = NULL;
  size_t max = 1024;
  XYZ *p = new XYZ[max]; 
  XYZ *p_Temp = NULL;   
  size_t nv = 0;
  int ntri = 0;
  double x, y;
  bool b_Ok = false;
    
  randomize();
  nv = 0;
  p = new XYZ[max];
  while (nv != n_MaxPoints){
    do{
      b_Ok = true;
	  x = (double) MISRENDER::frand() * 500.0;
	  y = (double) MISRENDER::frand() * 500.0;
      for(size_t n_Cpt = 0; n_Cpt <= nv; n_Cpt++){
        if((x == p[n_Cpt].x) && (y == p[n_Cpt].y)) b_Ok = false;
      }// to avoid similar points in the array
    }while(!b_Ok);
    if (nv >= max){
      max = max * 2;            // double the size of the array if necessary
      p_Temp = new XYZ[max]; 
      for (size_t i = 0; i < nv; i++) {
        p_Temp[i] = p[i];  
      }
      delete []p;  
      p = p_Temp; 
    }   
    p[nv].x = x * 1.0;
    p[nv].y = y * 1.0;
    nv++;
  }
  p_Temp = new XYZ[nv + 3]; 
  for (size_t i = 0; i < nv; i++) {
    p_Temp[i] = p[i];      
  }
  delete []p;           
  p = p_Temp;
  v = new ITRIANGLE[3 * nv]; 
  qsort(p, nv, sizeof(XYZ), XYZCompare);
  Triangulate(nv, p, v, ntri);
  
	std::vector<XYZ> pts;

	while( glfwGetWindowParam(GLFW_OPENED) )
	{
		glClear(GL_COLOR_BUFFER_BIT);
		outputtriangle(nv, p, v, ntri); // use this fonction to trace the mesh (via 
		
		glPointSize(2.0);
		glEnable(GL_POINT_SMOOTH);

		glBegin(GL_POINTS);
		for(size_t i=0; i<pts.size(); i++)
		{
			glVertex2d( pts[i].x, pts[i].y );
		}
		glEnd();


		std::stringstream ss;
		ss << "num pts: " << pts.size(); 
		glfwSetWindowTitle( ss.str().c_str() );

		glfwSwapBuffers();

		static bool canMake = true;

		if( pts.size() >= max )
		{
				canMake = false;
		}

		if( glfwGetMouseButton(GLFW_MOUSE_BUTTON_1) && canMake)
		{
			canMake = false;
			int mx,my;
			glfwGetMousePos(&mx,&my);
			XYZ newXYZ;
			newXYZ.x = mx;
			newXYZ.y = my;
			pts.push_back( newXYZ );
		}
		if( glfwGetMouseButton(GLFW_MOUSE_BUTTON_1) )
		{
		}else{
			canMake = true;
		}

		static bool canMake2 = true;
		if( glfwGetMouseButton(GLFW_MOUSE_BUTTON_2) && canMake2)
		{
			canMake2 = false;
			double r = 240;
			
			int numPts = 6*6+6;
			for(int circs=0; circs < 6; circs++)
			{
				r -= 30;
				numPts -= 3;
			
				for(int i=0; i<numPts; i++)
				{
					
					double ang = i * 2 * 3.14 / numPts;
					XYZ val;
					val.x = xres/2 + r * cos(ang) + MISRENDER::frand() * .1;
					val.y = yres/2 + r * sin(ang) + MISRENDER::frand() * .1;

					if( pts.size() < max )
						pts.push_back( val );
				
				}
			}
		}

		if( !glfwGetMouseButton(GLFW_MOUSE_BUTTON_2) )
		{
			canMake2 = true;
		}

		static bool canMake3 = true;
		if( glfwGetMouseButton(GLFW_MOUSE_BUTTON_3) && canMake3)
		{
			int wid = 25;
			int hei = 25;
			for(int i=0; i<wid; i++)
			{
				for (int j=0; j<hei; j++)
				{
					XYZ val;
					val.x = 50 + i * 17 + 0.001 * MISRENDER::frand();
					val.y = 50 + j * 17;

					if( pts.size() < max )
						pts.push_back( val );
				}
				
			}
		}

		if( !glfwGetMouseButton(GLFW_MOUSE_BUTTON_3) )
		{
			canMake3 = true;
		}


		if( glfwGetKey(GLFW_KEY_SPACE) )
		{
			if( pts.size() > 2 )
			{

			delete []p;
			delete []v;
			p = NULL;
			v = NULL;

			nv = 0;
			p = new XYZ[max];

			for(size_t i=0; i<pts.size(); i++)
			{
				p[nv].x = pts[i].x;
				p[nv].y = pts[i].y;
				nv++;
			}

			p_Temp = new XYZ[nv + 3]; 
			for (size_t i = 0; i < nv; i++) {
				p_Temp[i] = p[i];      
			}

			delete []p;           
			p = p_Temp;
			v = new ITRIANGLE[3 * nv]; 
			
			qsort(p, nv, sizeof(XYZ), XYZCompare);
			Triangulate(nv, p, v, ntri);

			pts.clear();
			}
		}

		
	}
	
  
	delete []p;// OpenGL, DirectX, ...)
  delete []v;
  p = NULL;
  v = NULL;


	glfwTerminate();
  return 0;
}




